 
C     $Id: cspline.F 17 2012-12-07 05:10:30Z wangsl2001@gmail.com $
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  2003  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     ##################################################################
c     ##                                                              ##
c     ##  subroutine isplpe  --  cubic periodic interpolating spline  ##
c     ##                                                              ##
c     ##################################################################
c
c
c     "isplpe" computes the coefficients for a cubic periodic
c     interpolating spline
c
c     literature reference:
c
c     G. Engeln-Mullges and F. Uhlig, Numerical Algorithms with Fortran,
c     Springer Verlag, 1996, Section 10.1.2
c
c
      subroutine isplpe (ng,xn,fn,mrep,b,c,d,h,du,dm,rc,rs)
      implicit none
      include 'iounit.i'
      integer i,ng,mrep,iflag
      real*8 eps,average
      real*8 temp1,temp2
      real*8 xn(0:ng),fn(0:ng)
      real*8 b(0:ng),c(0:ng)
      real*8 d(0:ng),h(0:ng)
      real*8 du(ng),dm(ng)
      real*8 rc(ng),rs(ng)
c
c
c     check the periodicity of fn, and for subsequent call
c
      eps = 0.00001d0
      if (abs(fn(ng)-fn(0)) .gt. eps) then
         write (iout,10)
   10    format (' ISPLPE  --  Warning, Input Values to Spline',
     &              ' are Not Periodic')
      end if
      average = 0.5d0 * (fn(0) + fn(ng))
      fn(0) = average
      fn(ng) = average
      if (mrep.ne.1 .and. mrep.ne.2)  return
c
c     get auxiliary variables and matrix elements on first call
c
      if (mrep .eq. 1) then
         do i = 0, ng-1
            h(i) = xn(i+1) - xn(i)
         end do
         h(ng) = h(0)
         do i = 1, ng-1
            du(i) = h(i)
         end do
         du(ng) = h(0)
         do i = 1, ng
            dm(i) = 2.0d0 * (h(i-1)+h(i))
         end do
      end if
c
c     compute the right hand side
c
      temp1 = (fn(1)-fn(0)) / h(0)
      do i = 1, ng-1, 1
         temp2 = (fn(i+1)-fn(i)) / h(i)
         rs(i)  = 3.0d0 * (temp2-temp1)
         temp1 = temp2
      end do
      rs(ng) = 3.0d0 * ((fn(1)-fn(0))/h(0)-temp1)
c
c     solve the linear system, with factorization on first call
c
      if (mrep .eq. 1) then
         call cytsy (ng,dm,du,rc,rs,c(1),iflag)
         if (iflag .ne. 1)  return
c
c     proceed without factorization on subsequent calls
c
      else
         call cytsys (ng,dm,du,rc,rs,c(1))
      end if
c
c     compute remaining spline coefficients
c
      c(0) = c(ng)
      do i = 0, ng-1
         b(i) = (fn(i+1)-fn(i))/h(i) - h(i)/3.0d0*(c(i+1)+2.0d0*c(i))
         d(i) = (c(i+1)-c(i)) / (3.0d0*h(i))
      end do
      b(ng) = (fn(1)-fn(ng))/h(ng) - h(ng)/3.0d0*(c(1)+2.0d0*c(ng))
      return
      end
c
c
c     #############################################################
c     ##                                                         ##
c     ##  subroutine cytsy  --  solve cyclic tridiagonal system  ##
c     ##                                                         ##
c     #############################################################
c
c
c     "cytsy" solves a system of linear equations for a cyclically
c     tridiagonal, symmetric, positive definite matrix
c
c     literature reference:
c
c     G. Engeln-Mullges and F. Uhlig, Numerical Algorithms with Fortran,
c     Springer Verlag, 1996, Section 4.11.2
c
c
      subroutine cytsy (ng,dm,du,cr,rs,x,iflag)
      implicit none
      integer ng,iflag
      real*8 dm(1:ng),du(1:ng)
      real*8 cr(1:ng),rs(1:ng)
      real*8 x(1:ng)
c
c
c     factorization of the input matrix
c
      iflag = -2
      if (ng .lt. 3)  return
      call cytsyp (ng,dm,du,cr,iflag)
c
c     update and backsubstitute as necessary
c
      if (iflag .eq. 1) then
         call cytsys (ng,dm,du,cr,rs,x)
      end if
      return
      end
c
c
c     #################################################################
c     ##                                                             ##
c     ##  subroutine cytsyp  --  tridiagonal Cholesky factorization  ##
c     ##                                                             ##
c     #################################################################
c
c
c     "cytsyp" finds the Cholesky factors of a cyclically tridiagonal
c     symmetric, positive definite matrix given by two vectors
c
c     literature reference:
c
c     G. Engeln-Mullges and F. Uhlig, Numerical Algorithms with Fortran,
c     Springer Verlag, 1996, Section 4.11.2
c
c
      subroutine cytsyp (ng,dm,du,cr,iflag)
      implicit none
      integer i,ng,iflag
      real*8 eps,row,d
      real*8 temp1,temp2
      real*8 dm(1:ng),du(1:ng)
      real*8 cr(1:ng)
c
c
c     set error bound and test for condition ng > 2
c
      eps = 0.00000001d0
      iflag = -2
      if (ng .lt. 3)  return
c
c     checking for positive definite matrix a and for strong
c     nonsingularity of a for ng=1 and if a is
c
      row = abs(dm(1)) + abs(du(1)) + abs(du(ng))
      if (row .eq. 0.0d0) then
         iflag = 0
         return
      end if
      d = 1.0d0 / row
      if (dm(1) .lt. 0.0d0) then
         iflag = -1
         return
      else if (abs(dm(1))*d .le. eps) then
         iflag = 0
         return
      end if
c
c     factoring a while checking for a positive definite and strong
c     nonsingular matrix a
c
      temp1 = du(1)
      du(1) = du(1) / dm(1)
      cr(1) = du(ng) / dm(1)
      do i = 2, ng-1
         row = abs(dm(i)) + abs(du(i)) + abs(temp1)
         if (row .eq. 0.0d0) then
            iflag = 0
            return
         end if
         d = 1.0d0 / row
         dm(i) = dm(i) - temp1*du(i-1)
         if (dm(i) .lt. 0.0d0) then
            iflag = -1
            return
         else if (abs(dm(i))*d .le. eps) then
            iflag = 0
            return
         end if
         if (i .lt. (ng-1)) then
            cr(i) = -temp1 * cr(i-1) / dm(i)
            temp1 = du(i)
            du(i) = du(i) / dm(i)
         else
            temp2 = du(i)
            du(i) = (du(i) - temp1*cr(i-1)) / dm(i)
         end if
      end do
      row = abs(du(ng)) + abs(dm(ng)) + abs(temp2)
      if (row .eq. 0.0d0) then
         iflag = 0
         return
      end if
      d = 1.0d0 / row
      dm(ng) = dm(ng) - dm(ng-1)*du(ng-1)*du(ng-1)
      temp1 = 0.0d0
      do i = 1, ng-2
         temp1 = temp1 + dm(i)*cr(i)*cr(i)
      end do
      dm(ng) = dm(ng) - temp1
      if (dm(ng) .lt. 0) then
         iflag = -1
         return
      else if (abs(dm(ng))*d .le. eps) then
         iflag = 0
         return
      end if
      iflag = 1
      return
      end
c
c
c     ################################################################
c     ##                                                            ##
c     ##  subroutine cytsys  --  tridiagonal solution from factors  ##
c     ##                                                            ##
c     ################################################################
c
c
c     "cytsys" solves a cyclically tridiagonal linear system
c     given the Cholesky factors
c
c     literature reference:
c
c     G. Engeln-Mullges and F. Uhlig, Numerical Algorithms with Fortran,
c     Springer Verlag, 1996, Section 4.11.2
c
c
      subroutine cytsys (ng,dm,du,cr,rs,x)
      implicit none
      integer i,ng
      real*8 sum,temp
      real*8 dm(1:ng),du(1:ng)
      real*8 cr(1:ng),rs(1:ng)
      real*8 x(1:ng)
c
c
c     updating phase
c
      temp = rs(1)
      rs(1) = temp / dm(1)
      sum = cr(1) * temp
      do i = 2, ng-1
         temp = rs(i) - du(i-1)*temp
         rs(i) = temp / dm(i)
         if (i .ne. (ng-1))  sum = sum + cr(i)*temp
      end do
      temp = rs(ng) - du(ng-1)*temp
      temp = temp - sum
      rs(ng) = temp / dm(ng)
c
c     backsubstitution phase
c
      x(ng) = rs(ng)
      x(ng-1) = rs(ng-1) - du(ng-1)*x(ng)
      do i = ng-2, 1, -1
         x(i) = rs(i) - du(i)*x(i+1) - cr(i)*x(ng)
      end do
      return
      end
